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We report a numerical calculation of the total number of disordered jammed configurations fl of 
N repulsive, three-dimensional spheres in a fixed volume V . To make these calculations tractable, 
we increase the computational efficiency of the approach of Xu et al. (Phys. Rev. Lett. 106, 245502 
(2011)) and Asenjo et al. (Phys. Rev. Lett. 112, 098002 (2014)) and we extend the method to 
allow computation of the configurational entropy as a function of pressure. The approach that we 
use computes the configurational entropy by sampling the absolute volume of basins of attraction 
of the stable packings in the potential energy landscape. We Hnd a surprisingly strong correlation 
between the pressure of a configuration and the volume of its basin of attraction in the potential 
energy landscape. This relation is well described by a power law. Our methodology to compute the 
number of minima in the potential energy landscape should be applicable to a wide range of other 
enumeration problems in statistical physics, string theory, cosmology and machine learning, that 
aim to find the distribution of the extrema of a scalar cost function that depends on many degrees 
of freedom. 


I. INTRODUCTION 

Many questions in physics are easy to pose but diffi¬ 
cult to answer. One such question is: how many micro¬ 
scopic states of a given system are compatible with its 
macroscopic properties? In statistical mechanics, knowl¬ 
edge of this number allows us to compute the entropy, 
and thereby predict the macroscopic properties of a sys¬ 
tem from knowledge of the interaction between atoms or 
molecules. 

In granular matter we can similarly ask how many mi¬ 
crostates are compatible with a given set of macroscopic 
properties. However, the computation of the correspond¬ 
ing absolute entropy has thus far proven to be extremely 
challenging. Without such knowledge, it is not possible 
to explore the analogies and differences between granular 
and Boltzmann entropy. Being able to compute the con¬ 
figurational entropy is therefore clearly important. The 
more so as granular materials are ubiquitous in every¬ 
day life (sand, soil, powders). Many industrial processes 
involve granular materials. In the natural world, the 
Earth’s surface contains vast granular assemblies such as 
dunes, which interact with wind, water, and vegetation 
[1] . Packings of particles that are soft or biological in na¬ 
ture, such as cells, hydrogels and foams are also known 
to undergo jamming [2] and their behaviour to be “gran¬ 
ular” viz not subject to thermal motion. Moreover, as 
glasses and granular materials share many properties it 
has been proposed that their physics may be controlled 
by the same underlying principles [3]. 

The study of granular materials is complicated by 
the fact that these materials are intrinsically out-of- 
equilibrium. In fact, thermal motion plays no role in 
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granular matter. It maintains its configuration unless 
driven by external forces. As a consequence, the proper¬ 
ties of granular materials depend upon their preparation 
protocol. 

Granular materials are athermal and cannot therefore 
be described by statistical mechanics. However, these 
materials can exists in a very large number of distinct 
states and this fact inspired Edwards and Oakeshott 
[4] well over two decades ago to propose a statistical- 
mechanics-like formalism to describe the properties of 
granular matter. In its original version, the Edwards 
theory assumed that all mechanically stable configura¬ 
tions (‘jammed’ states) are equally probable and that 
the logarithm of the number of these states plays a role 
similar to that of entropy. In this theoretical framework, 
the volume of the system and its compactivity (i.e. the 
derivative of volume with respect to the configurational 
entropy) are the analogues of the energy and temperature 
in thermal systems. 

In the absence of explicit calculations (or measure¬ 
ments) of the absolute configurational entropy, a direct 
test of the Edwards hypothesis has proven difficult, and 
different authors have arrived at different conclusions 
based on indirect tests in either simulations [5-8] or ex¬ 
periments [9, 10]. In addition, alternative definitions of 
entropy have been proposed to characterise the complex¬ 
ity of granular systems while circumventing explicit enu¬ 
meration of states [11, 12]. 

Numerous tests of the Edwards volume ensemble have 
focused on the determination of the compactivity [13-21]. 
However, the role of compactivity as a temperature-like 
quantity is problematic as Puckett and Daniels [22, 23] 
have shown that it does not satisfy the equivalent of the 
zero-th law of Thermodynamics - the law that is the basis 
of all thermometry. 

Edwards’ theory has been generalised to include the 
distribution of stresses within the system through the 
force-moment tensor [24-27] and another analogue of 
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temperature emerged, known as angoricity, which is a 
measure of the change in entropy with stress. The exper¬ 
iments by Puckett and Daniels [22] showed that angoric- 
ity, unlike compactivity, is a temperature-like quantity as 
it satishes the zero-th law. 

To date only a few examples of numerical tests of the 
generalised Edwards ensemble are available [22, 26, 28, 
29] . Numerical tests of the stress ensemble focus on sys¬ 
tems of soft spheres near jamming where the compactiv¬ 
ity X ^ oo and fluctuations in volume are negligible 
compared to stress fluctuations [26, 27]. Wang et al. 
[28, 29] proposed a unified test that compared ensemble 
averaged results over volume and stress with predictions 
for the jamming transition, finding agreement; we note, 
however, that in the latter approach the results rely sig- 
nihcantly on the equiprobability assumption. 

When the system is composed of very stiff grains, or 
is close to jamming, any small deformation will lead to 
a large change in the contact forces. In these limits the 
geometric and the force degrees of freedom can be de¬ 
coupled, giving rise to the force network ensemble [30] 
(FNE). In this framework, force networks are constructed 
on a given geometry and each force state is assumed to 
be equiprobable. The FNE has been utilised as a testing 
ground for statistical frameworks [31-33]. 

More than two decades after its introduction many fun¬ 
damental questions concerning the Edwards hypothesis 
remain unanswered. This unsatisfactory state of affairs 
is at least partly due to the fact that no efficient meth¬ 
ods existed to measure or compute the absolute conhg- 
urational entropy directly. Until recently, the only way 
to determine the conhgurational entropy was by direct 
enumeration of the distinct jammed states of a system. 
This method is inefficient and cannot be used for sys¬ 
tems that contain more than 10-20 particles. Over the 
past few years, the situation on the numerical front has 
changed: recent numerical work by Asenjo et al. [34, 35], 
based on an approach introduced by Xu et al. [36], has 
demonstrated that it is possible to compute the number 
of distinct jammed states of a system, even when this 
number is far too large (e.g. 10^®°) to allow direct enu¬ 
meration. The approach of Refs. [34-36] replaces an in¬ 
tractable enumeration problem by a tractable scheme to 
sample the (absolute) volume of the basins of attraction 
of stable states in the potential energy landscape. 

The approach described herein is completely general 
and it extends to any energy landscape problem that aims 
to find the extrema of a scalar cost function that depends 
on many degrees of freedom. Enumerating the number 
of solutions or stationary points, and their distribution, 
for certain classes of random functions is a classical prob¬ 
lem in mathematics and statistics [37-50]. In statistical 
physics, ad hoc numerical and theoretical methods have 
been developed in the realms of random Gaussian and 
polynomial fields [51-59]. In this sense, particular atten¬ 
tion has been devoted to the mean-field p-spin spherical 
model of a spin glass with quenched disorder [60-64]. 
A related area is the computation of the conhgurational 
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FIG. 1: Entropy as a function of system size N for 
two (Ref. [34]) and three-dimensional (this work) jammed 
sphere packings. Dashed curves are lines of best ht of the 
form S = aN. 


contribution to the entropy of structural glasses [65, 66]. 
The physical signihcance of this method goes even fur¬ 
ther to encompass string theory [67-69], cosmology [70- 
74] and machine learning [75-78]. 

We note that the geometrical structure of the basins of 
attraction of jammed states had been studied by O’Hern 
and co-workers [8, 79, 80]. O’Hern also reported direct 
enumeration estimates of the number of jammed states 
of small systems. A rather different technique (‘basin 
sampling’) to count the number of energy minima in the 
potential energy landscape of small clusters had been re¬ 
ported by Wales and co-workers [81, 82]. 

We note that, for the system (and protocol) considered 
by Asenjo et al., not all packings are equally probable. 
However, as shown in Ref. [34], the equal-probability hy¬ 
pothesis is not needed to arrive at a meaningful definition 
of an extensive granular entropy. When, in the remain¬ 
der of the present paper, we mention the configurational 
(granular) entropy, we refer to the definition of Ref. [34] . 

We stress that, even though the approach of Refs. [34- 
36] allows to solve enumeration problems which were far 
from possible using direct enumeration, it is still compu¬ 
tationally expensive. Thus far, it had only been applied 
to two-dimensional packings. Substantial ‘technical’ im¬ 
provements were needed to make the method fast enough 
to deal with three-dimensional systems. 

In the present paper, we present the first enumeration 
of the number of jammed packings for three-dimensional 
systems consisting of up to 128 soft spheres. A di¬ 
rect comparison of the entropy measured as a function 
of system size for two and three-dimensional jammed 
sphere packings is shown in Fig. 1. The potential of the 
method presented herein can be verified unequivocally 
from Fig. 1: we are able of tackling problems at least 
500 million times more complex, and of greater computa- 
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tional cost, than the already spectacularly difficult ques¬ 
tions confronted by Asenjo et al. [34]. Furthermore we 
show how our improved procedure allows first-principles 
computation of configurational entropy as a function of 
system size and pressure. The method and the technical 
improvements needed to overcome this numerical chal¬ 
lenge are presented alongside the main results. 

The remainder of this work is organised as follows. Sec¬ 
tion II describes the basic principle of the mean basin vol¬ 
ume method for counting, and explains how that strat¬ 
egy can be applied to enumerate granular packings. The 
enumeration and entropy results for three-dimensional 
jammed sphere packings as a function of system size 
and pressure are reported in Sec. III. The rest of the 
manuscript is dedicated to the description of the im¬ 
proved numerical method. Section IV outlines our pro¬ 
tocol for sampling different granular packings, and it de¬ 
scribes the corresponding potential energy landscape and 
minimisation techniques. Application of thermodynamic 
integration to compute the volume of a basin of attrac¬ 
tion in such a landscape is described in Sec. V. Aspects of 
the data analysis tools used on the histograms of sampled 
basin volumes, and related configurational entropy defi¬ 
nitions, are described in Sec. VI. Conclusions are drawn 
in Sec. VII. Further technical background is given in the 
appendices. 


II. BASIC PRINCIPLE: COUNTING BY 
SAMPLING 

In this section, we briefly review the numerical ap¬ 
proach that we use to compute the number of distinct 
jammed states. We stress that the approach that we use 
has much wider applicability than the counting of gran¬ 
ular packings [51-68, 70-78, 83]. In the context of gran¬ 
ular packings, our aim is to compute the number of ways 
O in which N spheres can be arranged in a given vol¬ 
ume Vbox of Euclidean dimension d. Knowledge of O al¬ 
lows us to compute configurational entropies and related 
quantities from first principles [4, 34]. Our approach is 
based on a rigorous mapping of the enumeration prob¬ 
lem onto counting the number of minima of a potential 
energy landscape [36]. The approach makes no use of a 
harmonic [84] or quasi-harmonic [85] approximation. For 
a system of hard particles the potential energy function is 
discontinuous, that is, the energy of the system is either 
zero, if no two particles overlap, or infinity otherwise. 
Then, at jamming, in the absence of rattlers, basins of 
attraction are single points in configuration space and 
they have no associated volume. This does not mean 
that we cannot sample the energy minima of a system 
of hard particles. The reason is that all jammed struc¬ 
tures of hard particles correspond to the zero potential 
energy minima of a system with a continuous repulsive 
potential with the same range as the hard-core diameter 
of the hard particles. In what follows, we focus on this 
class of systems, but we generalise the problem by also 


considering minima with a non-zero potential energy. In 
particular, we consider spherical particles with a hard 
core and a short-ranged continuous repulsive interaction. 
Under conditions where this system is jammed, a system 
with only the hard-core interactions would still be fluid 
and would sample the accessible configuration space uni¬ 
formly. This remaining accessible volume is partitioned 
in basins of attraction defined by the soft shells. The 
HS-WCA potential used to simulate hard-core plus soft- 
shell interactions and the packing preparation protocol 
are described in Sec IV B. For an illustration of the pack¬ 
ing preparation protocol refer to Fig. 2. As we argue 
below, using an HS-WCA model greatly improves the ef¬ 
ficiency of determination of basin volumes. 

Let us denote the total available volume in dN- 
dimensional space as V. Note that V is not the total 
volume of configuration space (U'^), but just that part 
of the volume that is free of hard-core overlaps. It is 
the configurational part of the partition function of the 
hard-core system at the number density under consider¬ 
ation. Since the accessible configuration space is tiled 
by the basins of attraction of the distinct energy minima 
[84, 86-88] we can write: 

n 

v = (1) 

where Vi is the volume of the i-th basin of attraction and 
n is the total number of distinct minima. We thus make 
the simple observation: 

a o ^ 

= H(u), (2) 

i—1 1 

where (v) is the mean basin volume, from which it follows 
immediately that 

V 

n = ^ (3) 

{v) 

We note that, for sphere packings, V is known from the 
equation of state of the underlying hard sphere fluid (see 
Appendix E) and we can measure (v) by thermodynamic 
integration, as discussed in detail in the Sec IV. The ap¬ 
proach of [34, 36] has thus turned the intractable enumer¬ 
ation problem of finding H into a sampling one, namely 
measuring {v). 

III. RESULTS: COUNTING DISORDERED 3D 

SPHERE PACKINGS 

The mean basin volume method for enumerating the 
number of mechanically stable packings was introduced 
by Xu et al. [36], and tested on a small system of soft 
disks. Asenjo et al. [34] then made a number of modifi¬ 
cations to the algorithm that allowed them to apply it to 
larger systems of up to 128 disks. As was the case with 
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FIG. 2: Hard sphere fluid at (/)hs = 0.5, left, and HS-WCA jammed packing at ^ss = 0.7, right, for a system 
of 44 polydisperse hard spheres with mean radius {vh) = 1 and standard deviation ctrs = 0.05. We prepare the 
polydisperse HS fluid configurations at fixed packing fraction (/)hs = 0.5 by a Monte Carlo simulation. Particles are 
then inflated by the same factor, proportional to their radius (spheres are coloured according to their radius), to obtain 
an over-compressed soft spheres jammed packing at ())gg = 0.7 by an infinitely fast quench (energy minimisation). 


Ref. [36], the calculations of Ref. [34] focused on two- 
dimensional systems because of the high computational 
costs involved in studying 3D systems. Here we present 
results for systems of three-dimensional soft spheres. We 
are thus in a position to compute the configurational en¬ 
tropy of a real (but idealised) three-dimensional system. 

We first describe an analysis similar to the one reported 
by Asenjo et al. [34] to verify the extensivity of the en¬ 
tropy S{V) at constant packing fraction. Next, we ex¬ 
tend our approach to the generalised Edwards ensemble, 
i.e. one based on a description of the system in terms 
of its volume and pressure, to compute the generalised 
entropy S{V,'P). 

We investigate three-dimensional packings with sys¬ 
tem sizes ranging from 24 to 128 HS-WCA particles, see 
Eq. (15), at (/iRS = 0.5 hard-sphere fluid packing frac¬ 
tion and ^gs = 0.7 soft sphere packing fraction, corre¬ 
sponding to a ratio of the soft and hard-sphere radii ratio 
rgg/rns = 1-12, prepared following the protocol outlined 
in Sec. IV. For each system size we compute the volume of 
the basin of attraction of approximately 1000 packings. 
Each PT run (see Sec. V) was performed on 15 paral¬ 
lel threads of a single eight-core dual Xeon E5 — 2670 
(2.6GHz, Westmere) node. Our choice of convergence 
criterion was such that when the uncorrelated statistical 
error for each of the replicas’ mean square displacement 
fell below 5% the calculations were terminated. This set¬ 
up translated in run times ranging from 10 to 300 hours 
per basin depending on system size, which amounts to 
0(10®) hours of total run time and O(IO^) total cpu 


hours. We then analyse the corresponding distribution 
of dimensionless free energies following the protocol de¬ 
scribed in Sec. V and VI and summarised in Appendix A. 


A. Extensivity of the entropy 

We first computed two alternative definitions of en¬ 
tropy: the Gibbs entropy Sq = — Pi lii(pi) — ln(V!) 
and Edwards (Boltzmann) entropy Sb = ln(D) — ln(iV!), 
where pi is the probability to sample packing i and D is 
the total number of mechanically stable states (or min¬ 
ima in the energy landscape). A detailed discussion of 
these definitions is outlined in Sec. VI. The results of 
these calculations are summarised in Fig. 3. Our results 
strongly suggest that, also in three dimensions, the en¬ 
tropy thus defined is extensive. Note that extensivity re¬ 
quires not only that the entropy scales linearly with sys¬ 
tem size, but also that it crosses zero at the origin. The 
slightly higher value of the Edwards entropy compared 
to the Gibbs entropy is consistent with the observation 
that Edwards’ equiprobability corresponds to the maxi¬ 
mum possible entropy of a system with D states. We also 
show that our estimates for the Edwards’ entropy are rel¬ 
atively insensitive to the precise strategy used to compute 
it. In Fig. 3, we compare three methods: a parametric fit 
to a generalised Gaussian cumulative distribution func¬ 
tion (c.d.f.) using a non-linear least squares method, a fit 
to the corresponding probability density function (p.d.f.) 
using maximum likelihood, and a non-parametric fit by 
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FIG. 3: Top left, entropy as a function of the system size 
N computed, in order, according to the Gibbs configura¬ 
tional entropy and the Edwards configurational entropy 
using a non-parametric fit by kernel density estimation 
(KDE), a parametric fit to a generalised Gaussian c.d.f. 
using a non-linear least squares method and a fit to the 
corresponding p.d.f. using maximum likelihood (ML). 
Gomparison of generalised Gaussian best-fit parameters 
for 2D (see Ref. [34]) and 3D sphere packings: scale pa¬ 
rameter tj (bottom left) and mean log-volume /i (top 
right) scale linearly with system size N; distributions are 
more peaked for 2D packings. In 2D we observe much 
stronger dependence of the shape parameter C (bottom 
right) as a function of system size than in 3D. 


kernel density estimation, which makes no a priori as¬ 
sumption about the shape of the distribution, other than 
the choice of the kernel function. We note, once again, 
that no post-processing is needed to compute the Gibbs 
version of the configurational entropy. Our results are in 
line with those reported by Asenjo et al. [34] for two- 
dimensional systems. 

The number of mechanically stable states O required 
by the Edwards’ definition of entropy is obtained subse¬ 
quently to fitting the numerically obtained distribution of 
log-basin volumes (dimensionless free energies) to a gen¬ 
eralised Gaussian distribution and unbiasing it appropri¬ 
ately, as described in Sec VI. We observe that the best-fit 
mean and scale parameters of the generalised Gaussian 
for the distribution of dimensionless free energies, ^ and 
a in Eq. (27) respectively, are also extensive, which al¬ 
though in line with what was found in two dimensions, 
is not a priori obvious. Finally we find that the shape 
parameter, C in Eq. (27), appears to depend only weakly 
on system size. The statistics are poor, but the data are 
compatible with the assumption that C 2 (Gaussian 
distribution) as N ^ oo. In 2D, the same limiting distri¬ 
bution of but with a much stronger size dependence, 
was observed. 


B. Entropy in the generalised Edwards ensemble 

We next consider the situation where the configura¬ 
tional entropy is a function of both the volume V and 
the stress tensor E of the system. The number of pack¬ 
ings with fixed V and E is denoted by (E, V). 

In the generalised Edwards ensemble [23, 26, 27], we fix 
the variables conjugate to V and E, viz. the compactivity 
X and the inverse angoricity tensor d. The generalised 
‘partition function’ can then be written as [23]: 

Edyn = (4) 

where Vv and are the volume and the force-moment 
(stress) tensor for state v. The weights w account for 
the protocol dependence of the probability to generate 
a state, and the sum runs over all mechanically stable 
states V. 

We can rewrite this partition function in terms of the 
density of states: 

^dy„= n // dS'"dVfIdyn(S,V^)e-'"/"'e-^(“^). 
i,k>r 

(5) 

Eor a system under hydrostatic pressure, and in the ab¬ 
sence of shear, we can write the force-moment tensor 
as E = /F, where F = VV = Tr(E)/3 is the internal 
Virial of the system. The inverse angoricity tensor d 
becomes a scalar a = dS/dV [27]. This result allows 
to simplify the notation significantly and at fixed vol¬ 
ume, through the mean basin volume method, we obtain 
the number of states integrated over all pressure states, 
n(V) = f dT’fl(V,T’). We now discuss how to gener¬ 
alise this procedure so that one can compute 
and therefore the configurational entropy, in the context 
of the generalised Edwards ensemble. 


1. Pressure to basin volume power-law relation 

To compute Z{X, a) directly, we would have to evalu¬ 
ate n(7^, V) as a function of both V and V. Whilst, with 
the tools that we have, this calculation is in principle 
possible, the computational costs would be several orders 
of magnitude larger than the, already quite substantial, 
costs of computing n(V). This would suggest that the 
computation of Z(X,a) is not possible at present. 

However, it turns out that we can still estimate the 
generalised configurational entropy because, as we dis¬ 
cuss below, we observe a surprisingly strong correlation 
between pressure and basin volume. 

From Fig. 4 we see that the basin volume for a given 
pressure state at fixed volume is strongly correlated with 
the pressure V. As the figure suggests, the relation be¬ 
tween — ln(u) and ln(7^) is approximately linear, and 
hence 

F(rix, </>ss) = - Hv) = - ln{r) + C{N), (6) 

K, 
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FIG. 4: Top: dimensionless free energy versus pressure 
of mechanically stable states at fixed volume for several 
system sizes. Best fit lines are in black. In the bottom 
left and right plots we show slope and intercept for each 
of the best fit lines as a function of system size. Both 
slope and intercept scale linearly with system size. 


where k, denotes the slope of the linear fit, and C{N) its 
value at 7^ = 1 (see Fig. 4). The value of k is not known 
a priori. It seems likely that k depends on the functional 
form of the potential. C{N) is a even less universal linear 
function of N, as it depends on the choice of units. 

We anticipate that this power law relationship survives 
for packings in two dimensions for a wide spectrum of 
packing fractions 4> > 4>j viz. as long as the system is 
jammed and sufficiently over-compressed [89]. 


2. Gibbs configurational entropy 

Using our approximate relation between pressure and 
basin volume, we can now rewrite Eq. (6) in terms of the 
probabilities for each jammed state 

ln(p,) = --ln(U.)-C(fV)-ln(V), (7) 

K 

which when substituted in the definition for the Gibbs 
entropy Eq. (22), gives the configurational entropy at a 
given volume in terms of the biased mean log-pressure 

Sg = - (ln(U))g + C{N) + ln(V) - ln(iV!). (8) 

K 

The significance of this equation should be apparent: for 
a sufficiently over-compressed packings of soft spheres at 



FIG. 5: Empirical cumulative distribution functions of 
the pressures for several system sizes. Dashed lines in 
the corresponding colour are curves of best fit to a gen¬ 
eralised log-normal distribution. The curves are mostly 
indistinguishable. Inset: best fit parameters for the gen¬ 
eralised log-normal distribution as a function of system 
size. The mean /r and scale parameter a scale linearly 
with 1/N, while the shape parameter C is approximately 
insensitive with respect to system size. 


a given packing fraction, the Gibbs configurational en¬ 
tropy can be approximately computed from sole knowl¬ 
edge of the average pressure, provided that n is known. 


3. Generalised Edwards configurational entropy 

To recover the number of states as a function of volume 
and stress we note that 

n{v,v) = n{v) u{x\v)dx, (9) 
Jv 

where U{V\V) is the unbiased probability distribution 
function of stresses at some specified volume. The di¬ 
rectly measured distribution of pressures depends on the 
protocol with which packings are generated. 

We distinguish between the biased, B{V\N,(j)ss) (as 
sampled by the packing protocol), and the unbiased, 
U{V\N, (j)ss), pressure distributions. Since the configu¬ 
rations were sampled proportional to the volume of their 
basin of attraction, using Eq. (6) we can compute the 
unbiased distribution analogously to Eq. (25) as 

U{r\N,fiss) = Q{N,fiss)B{r\N,cl)ss)e^^^'>V^^^, (10) 

where Q{N, (j)ss) = (j)ss) is the normalisation con¬ 

stant. 

Upon substitution of lnr2(U) = ln(V) — ln((u)) and 
of Eq. (10) for U{x\V), we write an expression for the 
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FIG. 6: Top: generalised Edwards entropy at fixed 
volume fraction for various system sizes. The curves 
show a well defined maximum for all sizes, while their 
shape depends on the specific parameters of the gener¬ 
alised log-normal that best fits the underlying distribu¬ 
tion of pressures. Bottom left: comparison between the 
Edwards entropy and the maximum value attained by 
each curve: max[5'B(7^, E)] scales linearly with size and 
its value is progressively closer to the marginal (total) 
Edwards entropy S'b(E), consistent with the fact that 
SsiV, V) is a negative exponential function, and the area 
under the curve is dominated by the mode for increasing 
system size. Sb{V) should constitute an upper limit to 
max[5'B(7^, E)] < Sb{V) and the two should be equiv¬ 
alent only in the thermodynamic limit. Bottom right: 
ensemble average of the pressure computed as a function 
of inverse angoricity a and system size. The curves, in 
the same colour as the top figure, do not diverge and the 
arrows indicates their value at a = 0. 


Edwards entropy as a function of volume and pressure 

/ i^v+sv 

SB{V,r) =\nlj Bix\V)x^/'^dx 

+ \n{V)+C{N)-ln{m). 

We fit the empirical cumulative distribution function 
(c.d.f.) of B{V) with the generalised log-normal c.d.f. cor¬ 
responding to Eq. (28) (see Fig. 5). We then numerically 
evaluate the generalised Edwards entropy SBi'PjV) at 
fixed volume, as shown in Fig. 6. 

In the thermodynamic limit we find 

Ssi'Pss) = 1 + c + -—- ln((/)ss) - /ex(<(>Hs), (12) 



where c = C{N)/N, see Appendix F for details of the 
derivation and further discussion. 

In Fig. 6 we also show the predicted expectation value 
for the pressure obtained via the ensemble average at 
arbitrary inverse angoricity a, 


(V) 


(ens) _ J 0 


VB{V\V)V^/^e-°‘^^ dV 


BiV\V)V^^'^e-°‘^^ dV 


(13) 


IV. PACKING PREPARATION PROTOCOL 
A. Sampling packings 


The physical properties of granular packings may de¬ 
pend strongly on the preparation protocol. This is il¬ 
lustrated by the Lubachevsky-Stillinger algorithm (LSA) 
procedure to prepare jammed packings of hard parti¬ 
cles [90] by compression (or, equivalently, by ‘inflation’ of 
the particles). If a monodisperse HS fluid is compressed 
rapidly the LSA will generate a low volume-fraction dis¬ 
ordered packing. However, for (very) slow compression 
rates, LSA will produce dense crystals [90, 91]. 

In the present work, we study a fluid of polydisperse 
spheres. We use a protocol related to a Stillinger-Weber 
quench that maps each fluid state to a local minimum, 
or “inherent structure”, connected by a path of steepest 
descent [86, 92]. 

To prepare the polydisperse fluid, we draw N 
particle radii {r}Ar from a Gaussian distribution 
Normal(l, cths) > 0, truncated at r = 0 (note that in 
our application the standard deviation cths is sufficiently 
small that it is extremely improbable to ever sample a 
negative radius). We set the box size to meet the target 
packing fraction of the hard sphere fluid ^hs and then 
place the particles in a random valid initial hard spheres 
configuration. The initial configuration is then evolved 
by a MC simulation [93] consisting of single particle ran¬ 
dom displacements and particle-particle swaps, and after 
equilibration, new configurations are recorded at regular 
intervals. We choose the length of these intervals such 
that, on average, each particle diffuses over a distance 
equal to the diameter of the largest particle. As long as 
(fins is well below the volume fraction where the fluid un¬ 
dergoes structural arrest, the allowed configurations of 
the fluid can be sampled uniformly. Importantly, this 
volume fraction is well below the random close packing 
« 0.64 and « 0.82 [79]). 

Given these HS fluid configurations, we now switch 
on the soft, repulsive interaction to generate over¬ 
compressed jammed packings of the particles (see Fig. 2). 
The particles are inflated with a WCA-like potential [94] 
to reach the target soft packing fraction (f>ss > > 

(fns- The hard spheres are inflated proportional to their 
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radius, so that the soft sphere radius is 


rs = 



\jd 

rh, 


(14) 


where d is the dimensionality of the box. Clearly, this 
procedure does not change the polydispersity of the sam¬ 
ple. 


B. Soft shells and minimisation 


We define the WCA-like potential around a hard core 
as follows: consider two spherical particles with hard core 
distance rh and soft core contact distance = rh{l + 
0), with 0 = (((>ss/</>Hs)^^‘' - 1- We can then write a 
horizontally shifted hard-sphere plus WCA (HS-WCA) 
potential as 


whs-wca(?’) 


' 00 
4e 

< 


0 



r < Th, 


Th <r < rs, 


r >rs 


(15) 

where a{rh) = {20 + 0'^)r\l2}l^ guarantees that the po¬ 
tential goes to zero at r^. For computational convenience 
(avoidance of square-root evaluations), the potential in 
Eq. 15 differs from the WCA form in that the inter¬ 
particle distance in the denominator of the WCA po¬ 
tential has been replaced with a difference of squares. 
Note that this implies that our potential resembles a 6-3 
potential more than a 12-6 potential. For our purpose, 
this difference is immaterial: we just need a short-ranged 
repulsive potential that diverges at the hard-core diame¬ 
ter and vanishes continuously at the soft-core diameter. 
The functional form of this potential is very similar to 
the HS-WCA potential used by Asenjo et al. [34], but 
cheaper to compute. We note that this potential is a 
type function, that is, its first derivative is continuous 
but not differentiable and its second derivative is discon¬ 
tinuous at Tfi. We take advantage of this property for the 
identification of rattlers (non-jammed particles) in our 
packings. 

Numerically evaluating this potential, we match the 
gradient and linearly continue the function uhs-wca(?') 
for r < r/j -|- £, with £ > 0 an arbitrary small constant, 
such that minimisation is still meaningful if overlaps do 
occur. 

The HS-WCA pair-potential was implemented using 
cell-lists [95, 96] with periodic boundary conditions, guar¬ 
anteeing 0{N) time complexity to the energy and gra¬ 
dient evaluations. Energy minimisations were performed 
with the CG_DESCENT algorithm [97-99] which, com¬ 
pared to FIRE [96, 100], reduces the average number of 
function evaluations for our system by a factor of 5, while 
preserving many of its desirable properties. 


The basin of attraction of a given minimum-energy 
configuration is the collection of all points connected to 
that minimum via a path of steepest descent [81, 101]. 
To measure the volume of a basin of attraction in the 
PES, we use thermodynamic integration [102, 103] and 
parallel tempering (PT) [93, 104-106]. 

The basic idea behind the method is that, but for the 
sign, the logarithm of the basin volume can be viewed as 
a dimensionless free energy. We cannot determine this 
free energy directly. We now switch on an increasingly 
harmonic potential that has its minimum at the mini¬ 
mum of the basin. In the limit of very large coupling 
constants (how large depends on the shape of the basin) 
the boundaries of the basin no longer affect the free en¬ 
ergy of the system, which has effectively been reduced 
to a dN dimensional harmonic oscillator with known free 
energy (for more details, see Appendix D). For zero cou¬ 
pling constant, instead, the system is completely uncon¬ 
strained and therefore in the state of interest. Thermo¬ 
dynamic integration allows us to compute the free energy 
difference between a reference state of known free energy 
and the (unknown) free energy associated with the orig¬ 
inal basin of attraction. 

A closely related approach is often used to compute the 
free energy of crystals of particles with a discontinuous 
potential, such as hard spheres [102, 103, 107]. Details of 
that method are summarised in the Appendix B, and the 
extension of the technique to basin volume measurement 
is described below. Details of the Hamiltonian PT are 
discussed in Appendix C. 


A. Free energy calculation for basin volumes 

To measure the volume of a basin by thermodynamic 
integration, we perform a walk inside the basin, that is, 
we start the MCMC random walk from the minimum 
energy configuration and we reject every move that 
takes us outside the basin [34, 36, 83]. This procedure 
can be cast in normal Monte Carlo language by defining 
an effective potential energy function (oracle) [/^(rjri) 
which is zero inside the basin and infinite outside. We 
can then write the volume of the basin: 


Vi = 




(16) 


In order for the oracle to test whether a proposed config¬ 
uration is inside or outside the basin, a full energy min¬ 
imisation must be performed. The numerous potential 
energy calls required for a full energy minimisation rep¬ 
resent the major obstacle to the scalability of the method. 

We view the negative log-basin-volume as a dimension¬ 
less free energy F( = — In(ui) [36] and compute it by 
thermodynamic integration, as described in Appendix B. 








9 


Therefore, we write, analogously to Eq. (B2): 

r^max 

In = Fhari^mayi) 2 J 

where denotes the coordinates of the i-th energy mini¬ 
mum. Unless k^ax, the maximum spring constant of the 
harmonic reference system, is very large, a finite fraction 
of the points belonging to the purely harmonic reference 
system will be located in the region where Ub = oo. 

We can correct for this effect in our calculation of 
Ehar(fcmax) by Computing the ratio of the partition func¬ 
tions of a system with a harmonic spring constant fcmax, 
both with and without the basin potential energy func¬ 
tion Ub- This ratio is given by 


7^ = 


y drexp[-U(r|ri,fcinax) - UB{r\ri)] 
/ drexp[-U(r|ri,fcmax)] 


(18) 


where V is the sum of the hard-core potential and 
the harmonic potential with spring constant fcmax, see 
Eq. (Bl). We note that TZ can be computed using a 
‘static’ (i.e. non-Markov chain) Monte Carlo simulation, 
sampling directly from the Boltzmann distribution of the 
harmonic oscillator with spring constant fcmax- Since the 
integral in the denominator is known [see Eq. (D2)], we 
write the dimensionless free energy of the harmonic ref¬ 
erence state for basin i as 

Fharikmax) = “^ In ) “InT^. (19) 

^ \ ^max / 


We note that, in order to avoid a singularity in the 
integrand, it is useful to perform the simulations fixing 
the centre of mass. It follows that the same corrections 
to the free energy as derived in Refs. [102, 103, 107] must 
be applied: similarly to Eq. (B5), but with the additional 
correction in Eq. (19), we write the basin volume as: 


-In?;, = AF(™)-ln(Ubox) 

(iV- l)d, / 27r 
-In ' 


— In 7^, 


( 20 ) 


where is the integral in Eq. (17), and the en¬ 

semble averages have been computed with a constrained 
centre of mass and it is evaluated as in Eq. (B13). 

Eigure 7 shows an example of the mean squared 
displacement (|r —ro|^)^, as a function of the spring 
constant k, along with the approximate expression in 
Eq. (B8) used to construct the change of variables in 
Eq. (B13). The resulting integrand, after the variable 
transform, is shown in the inset of Fig. 7. 


VI. BASIN VOLUME DISTRIBUTIONS AND 
DATA ANALYSIS 

Once the volumes of multiple basins have been sam¬ 
pled, these data can be used to compute the number of 


a 





- o - ■;o----o-n 


0 100 200 300 400 500 600 700 800 
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FIG. 7: Average squared displacement (|r — ro|^)^asa 
function of the spring constant k (symbols). The dashed 
line shows the expression in Eq. (B8). The data is mea¬ 
sured for a packing of A = 32 spheres, with (j)BS = 0.5 
and (pss = 0.7 via Hamiltonian PT. Inset: corresponding 
integrand for the thermodynamic integration, resulting 
from the change of variables in Eq. (B13). 


distinct packings [34], and from that, the Edwards en¬ 
tropy [4]. Furthermore we analyse the distribution of 
pressures of the different energy minima at given volume. 
In this work, we express pressure and volume in reduced 
units V/V* and v/v* everywhere with v* = (47r/3)(r^) 
and V* = ejv* being the units of volume and pressure, 
respectively. 


A. Gibbs configurational entropy 

Let us first consider the ‘Gibbs’ configurational en¬ 
tropy, Sg, defined by Asenjo et al. [34]: 

n 

Fg =-^Piln(p,)-ln(A!), (21) 

1=1 

where pi is the probability to sample packing i. For our 
preparation protocol, packings are sampled according to 
the volume of their basin of attraction, such that Pi = 
Vi/V. Then Eq. (21) gives 

n 

5G = -^p,ln(?;,)-kln(V)-ln(A!) 

= (F)B-hln(V)-ln(A!). 

The sum in Eq. (22) is the mean of the negative log- 
basin volumes (dimensionless free energies), as computed 
above, and weighted by the probabilities of preparing 
the corresponding basins. Therefore, the entropy can be 
obtained directly, and without approximation, from the 
sampled mean basin dimensionless free energy. 

















10 


From Eq. (22) we can also write the entropy per par¬ 
ticle in the thermodynamic limit as 

Sg(<('Ss) = 1 + (/)b + ln(0ss) - /ex(^Hs), (23) 

where /ex(<('Hs) is the excess free energy of the hard 
spheres fluid. In deriving this results we used Stirling’s 
approximation for large N and the fact that Vbox/v* = 

N/(t)ss- 


B. Edwards configurational entropy 

Edwards [4] suggested a Boltzmann-like entropy, where 
S equals the logarithm of the total number of pack¬ 
ings. Asenjo et al. [34] showed that, even for poly- 
disperse particles, indistinguishability of macrostates re¬ 
quires that 


S'b = In(fl)-ln(Af!), (24) 

The subtraction of ln(A^!) is necessary to guarantee ex- 
tensivity of the entropy. Unlike the Gibbs definition 
of entropy, Eq. (24) makes the explicit assumption of 
equiprobability of states. 

Eor a direct computation of the number of packings 
n, using Eq. (3), we need the average basin volume 
{v). Since our preparation protocol samples each mini¬ 
mum with a probability proportional to the volume of its 
basin of attraction, our samples of v are biased accord¬ 
ingly. Therefore, to obtain the unbiased average basin 
volume (w), the sampled basin volume distribution needs 
to be unbiased [34, 35, 83]. The unbiasing method used 
in the following work requires an analytical (or at least 
numerically integrable) description of the biased basin 
free energy distribution function. Different approaches 
to modelling this distribution give rise to somewhat dif¬ 
ferent analysis methods, which all yield consistent results. 
Again, we stress that no such additional steps are needed 
to compute the ‘Gibbs’ version of the configurational en¬ 
tropy. 

We distinguish between the biased, B{F\N,(j)ss) (as 
sampled by the packing protocol), and the unbiased, 
U{F\N,(j)^s)^ free energy distributions. Since the con¬ 
figurations were sampled proportional to the volume of 
their basin of attraction, we can compute the unbiased 
distribution as 


U{F\N, ^ss) = Q{N, ^ss)B{F\N, ^ss)e^ (25) 


where Q{N,(j)ss) is the normalisation constant 


QiN,^ss) 


dFBiF\N,(l>ss)e^ 




{v){N, (j)ss)- 


(26) 

Erom Eq. (25), unbiasing of the raw free energy distri¬ 
bution seems straightforward, however Asenjo at al. [34] 
noted that the most probable basins are about O(IO^) 
more probable than the small ones. Upon unbiasing, this 


factor is multiplied by a factor of about hence they 
observe that small basins are much more numerous than 
large ones and grossly under-sampled. 

To overcome this problem, one can fit the biased mea¬ 
sured free energy distribution B{F\N, (pss) and perform 
the unbiasing via Eq. (26) on the best fitting distribu¬ 
tion. B{F\N,(j)ss) must be bounded, hence it should de¬ 
cay with a functional form exp(—E*^) where v > 1. 

Before performing the fit we remove outliers from the 
free energy distribution following the distance-based out¬ 
lier removal method introduced by Knorr and Ng [108]. 
This is a form of clustering for which we choose to keep 
only those points for which at least half of the remain¬ 
ing data set is within 3 ct from the point, where a is the 
standard deviation computed for the raw data set. This 
procedure typically results in the exclusion of one or two 
points and it is essential for a successful fit to a gener¬ 
alised Gaussian model. 


1. Generalised Gaussian 


Assuming that U{F\N, (pss) is unimodal, which has 
been verified for very small systems [36], one can fit 
the raw distribution B{F\N, (pss) with a three-parameter 
generalised normal distribution 


PiF\F,a, C) 


2ar(l/C) 


(^) 


(27) 


where r(x) is the gamma function, a is the scale parame¬ 
ter, C is the shape parameter and F is the mean (free en¬ 
ergy) with variance tT^r(3/C)/r(l/C). In the limit C 2 
we recover the Gaussian distribution with standard devi¬ 
ation a. In practice it appears to be most stable to fit the 
empirical biased cumulative distribution function, rather 
than the histogram shape [34]. Alternatively, we also 
tested fitting to the observed p.d.f. with the maximum- 
likelihood method, obtaining consistent, but more scat¬ 
tered, results (see also Sec. III). 


2. Kernel density estimate 

To relax the assumption that the empirical distribu¬ 
tions can be fitted by a generalised Gaussian, one can 
also describe the distributions by kernel density estima¬ 
tion [109, 110]. Bandwidth selection is then done using 
Silverman’s rule of thumb as the initial guess for inte¬ 
grated squared error cross-validation [111]. The numer¬ 
ical integration step is performed, as for the generalised 
Gaussian description, via Eq. (26). 


C. Distribution of pressures 

In Sec. IIIB we have established a link between the 
pressure of a packing and the volume of its basin of at- 
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traction. In order to compute the entropy as a function 
of volume and pressure it is necessary to unbias the dis¬ 
tribution of pressures with respect to the sampling bias 
exp(—F), analogous to the previous section. We choose 
to describe the distribution of pressures V using the gen¬ 
eralised log-normal distribution [112] 


p(F’|ln(F),cr, C) 


2(C-Hi)/Car(l/C) 



ln(F) - ln(F) 

(7 



(28) 


with the first term on the r.h.s. being the normalisation 
constant and the remaining notation analogous to that 
of Eq. (27). For C = 2 this distribution reduces to the 
log-normal distribution. 


VII. CONCLUSIONS 


see for instance the discussion of the stress ensemble in 
the limit (j) ^ (j)j hy Henkes and Chakraborty [26, 27] 
and the work on the force network ensemble for systems 
of almost hard grains [31-33]. 
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The study of a statistical mechanics of granular materi¬ 
als has been complicated by the impossibility of directly 
computing fundamental thermodynamic quantities. In 
the present paper we have shown that configurational 
entropies of three-dimensional packings can, in fact, be 
computed. 

We have presented a method for the direct enumera¬ 
tion of the mechanically stable states of systems consist¬ 
ing of up to 128 frictionless soft three-dimensional spheres 
and we have shown that a definition of extensive entropy 
is possible, in line with the results for two dimensional 
systems reported by Asenjo et al. [34], with very minor 
differences in our observations. The study of 3D pack¬ 
ings is computationally demanding: the computational 
time required for each packing ranged between 10 and 
10^ cpu hours, depending on system size. The present 
study therefore required substantial algorithmic optimi¬ 
sation. 

We find that there is an approximately linear relation¬ 
ship between the logarithm of the pressure of a mechan¬ 
ically stable configuration and the logarithm of the vol¬ 
ume of its basin of attraction. 

The unexpected power law relationship between pres¬ 
sure and basin volume provides a way to extend our ap¬ 
proach to the generalised Edwards ensemble. We can 
analytically unbias the observed distribution of pressures 
and compute the entropy as a function of pressure at 
a given volume. Hence we have obtained consistent ex¬ 
pressions for the entropy in the thermodynamic limit. 
Knowledge of this distribution enables the first direct 
computation of angoricity. 

Tackling the study of granular materials from the en¬ 
ergy landscapes point of view is rather advantageous, al¬ 
though this does not come without burdens. This sort 
of approach is limited to soft frictionless particles, and 
we expect it to be reliable only at 4> > (j)j when the 
system is at least slightly over-compressed. Other theo¬ 
retical approaches are useful in more limiting situations. 


Appendix A: Basin volume method summary 

In summary, to count the number of ways spheres can 
pack into a given volume, we use the mean basin volume 
method outlined in Sec. H. We perform the following sim¬ 
ulations and analysis steps to obtain the required results: 

1. Obtain a number of different snapshots of an equili¬ 
brated hard sphere fluid at the desired volume frac¬ 
tion as described in Sec. IV A. This procedure 
fixes the number of measured basin volumes. 

2. Over-compress the sphere configuration by adding 
a soft shell. This compression yields, upon energy 
minimisation, a jammed packing with soft volume 
fraction ())ss > 

3. Estimate the maximum spring constant for the 
PT simulations, fcmax in Eq. (17), such that p in 
Eq. (18) reaches a value between 0.85 and 0.9. 
This is done by direct sampling and also gives the 
value of the average squared displacement for fcmax, 
(k-rol^) . 

'^max 

4. Obtain a preliminary estimate of the average 
squared displacement without harmonic tethering, 
(jr — rol^)g, by performing a MGMG walk in the 
basin. Use this result, with the estimate of /cmax 
from the previous step, to determine the spring con¬ 
stants k for the PT simulation, using Eqs. (Bll) 
and (B14). 

5. Perform a PT simulation to sample (jr — rg]^)^, as 
described in Appendix G. 

6. Compute the volume in Eq. (20) for each basin and 
analyse the distributions for all basins, at fixed vol¬ 
ume fraction and number of particles, as discussed 
in Sec. VI. This makes use of the total accessible 
volume, computed in Appendix E. 
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Section III shows examples of the type of results that can 
be obtained. Evaluation and minimisation of potential 
energy functions was performed with the pele [96] and 
PyCG_DESCENT [99] software packages. Monte Carlo 
simulations were performed with the mcpele package [93] . 


Appendix B: Free energy calculation for solids 

To compute the free energy of a system with discontin¬ 
uous potential energy function (e.g., hard disks or hard 
spheres), we construct a reversible path to the corre¬ 
sponding Einstein solid (see e.g. [107]). The harmonic 
potential with spring constant k is switched on while 
maintaining the hard core interactions intact: 

E(rlro,fc) = VHs(r) + A:Ehar(i’|ro) 

= VHs(r) + ^fc|I■-rol^ ^ ^ 

where Tq are the equilibrium coordinates of the Einstein 
crystal and VHs(r) denotes the hard core interactions. 
We can then compute the free energy difference between 
the Einstein crystal and the hard core system by evalu¬ 
ating the integral: 


pRS — Phar(^max) 



gE(rlro,fc) \ 

9 k Ik 


(B2) 


As discussed in Appendix D, we take the centre of mass 
to be fixed to avoid numerical issues in the limit fc —> 0. 
For a system with fixed centre of mass, we write the free 
energy difference between the target and the reference 
state as 


^^(CM) ^ ^(CM) _ 

From the partition function of the Einstein crystal with 
fixed centre of mass, Eq. (D7), and for the unconstrained 
crystal, Eq. (Dll), we can rewrite Eq. (B3) and rearrange 
it for the free energy of the unconstrained crystal: 


F = + In (iP(rcM = 0)) 


In 


/27r^^ 

V kiYisLK 



(B4) 


where the last term is Ekar and the second and third 
terms on the right hand side are the CM corrections for 
the unconstrained and the constrained solid, respectively. 
For a system with unit cell identical to the simulation box 
(with periodic boundary conditions), we have V{rcM = 
0) = l/14ox- Assuming that all particles have unit mass 
we can rewrite Eq. (B4) as 

F = AF^^^^ - In (Ebox) - In (^) ■ (B5) 

We are only left with AF^'^^I which can be found by 
evaluating the integral in Eq. (B2). In order to do so, we 


would like the integrand to be a well behaved function, 
possibly flat, permitting Gauss-Lobatto (CL) quadrature 
[113]. We transform the integration variable so that 




^ r\h 1 

(^rnax) 


2\CCM) 
k 


/G-^O) 


d[G-\fc)] g(fc)l(|r-ronf“^ 


(B6) 

where gik) is some function of k and G ^{k) is the prim¬ 
itive of the function 1/g(k). 

To choose an appropriate g{k), we note that in 
Eq. (D8) for very large k the mean squared displacement 
for the solid is 

(iV- l)d 


(Ir-i-ol = 


(B7) 


For k other than fcmax, we expect the mean squared dis¬ 
placement to depend on some effective spring constant. 
Hence we write 


(|r-ro|\ 
ed 


{N-l)d 


(B8) 


(B9) 


(k + C) 

such that the mean squared displacement at A: = 0 is 

{N-l)d 

e ’ 

from which we find ^ = {N — l)ci/(jr — rol^)fe=o [note 
that we can self consistently replace this definition for ^ 
in Eq. (B8) to obtain an approximation for the mean 
squared displacement at arbitrary k]. We would like 
the integrand g{k){\r — ro\‘^)k in Eq. (B6) to be roughly 
constant. Given the considerations above we choose 
g{k) Ki k + ^. One can easily verify that the integrand 
is now approximately constant. We can then rewrite the 
integral in Eq. (B6) as 

/■In(femax-I-C) 

^^(CM) ^ / 


/ln({) 




d [ln(A: -|- ^)] 


(BIO) 

Finally, to integrate Eq. (BIO) by GL quadrature, we 
require a variable, t, such that the integral upper and 
lower bounds are [—1,1]: 

21n(l-hfc/C)-l 


t = 


with differential 
dt = 


ln(l -bfcmax/C) 

d [In (1 -h /c/C)] . 


ln(l -h /Cmax/C) 

Therefore we rewrite Eq. (BIO) as a function of t: 

k. 


(Bll) 


(B12) 




dt In 1-1- 


/-il V C 

[fc(t)+C]^(|r-ro|2)f^^ 


(B13) 
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where k{t) can be found by rearranging Eq. (Bll). An 
example of the variable transform is shown in Fig. 7. 

It is straightforward to perform GL quadrature for a 
general number of abscissas n > 2 [113], because 

1 n— 1 

dtf{t) = Wif{-1) + + Wnf{l), (B14) 

-1 . — 9 


where rg denotes the equilibrium position, the index i 
denotes the Ath of N particles, each with d degrees of 
freedom, and we have assumed that the spring constant k 
is the same for all directions of motion. We can compute 
the mean squared particle displacement for a harmonic 
oscillator in the canonical ensemble analytically. We start 
with the partition function: 


where Wi are the weights and ti are the abscissas. The 
abscissas different from —1,1 are the n — 2 roots of 
dPn-iit)/dt, with Pn-i a Legendre polynomial. We 
evaluate this sum numerically using Numpy’s Legendre 
module [114]. The weights Wi can also be evaluated nu¬ 
merically for general n > 2, since they are related to 
Pn-i evaluated at ti [113]. For all results in this work, 
we choose n = 16 abscissas. 


Appendix C: Sampling the integrand: Hamiltonian 
Parallel Tempering 

To compute the integral in Eq. (B13), we need to mea¬ 
sure the integrand for different values oi k, as given by 
Eq. (Bll). Equilibration of the corresponding simula¬ 
tions can be accelerated using extensions of the parallel 
tempering technique, where replicas differ in chemical po¬ 
tential [115] or in the potential energy function [116, 117]. 

The Parallel-tempering acceptance rule for a swap of 
configurations between replicas with different Hamiltoni¬ 
ans follows from the condition of detailed balance: 


= {%) • 

We consider the free energy for the system F = 
InZ and observe that 


fdF{k) 

dk 


NVT 


f) 37 

O -| 

2 


/ 

J —c 


^-I3k\r-raf /2 


1 

2 


(|r 


roP), 


(D3) 

hence we can compute the mean squared distance for a 
dA-dimensional harmonic oscillator 


(|r-ron=2 


fdF{k) 
\ dk 


NVT 


dN 

Jk' 


(D4) 


acc[(r,, V,), irj,Vj) -» (rj, E,), (r,, Vj)] 
acc[(r„ Vj), {rj,Vi) (r„ V,), (r^, Vj)] 

^ exp{-^[V)(rj) -tEj-(ri)]} 

exp{-/3[V,{r,) + Vj{rj)]} 

= exp {-(3 [(E,(rj) -t V,-(ri)) - (Ei(r,) -t V,-(rj))]} , 

(Cl) 

where acc[- —?► •] denotes the swap acceptance probability. 
For the particular case of replicas coupled to a reference 
state rg by a harmonic potential with different coupling 
strengths ki, we find the swap acceptance rule 

acc[(r,, Vi), irj,Vj) (r,, V,)] 

= min |l, exp[^ [{kj - h){\rj - rgj^ - [r, - rgj^)] | . 

(C2) 

To check whether the replicas are well equilibrated, we 
consider the correlations in the “time series” of jr — rgjfc 
vs number of Monte Carlo steps for each replica. 


Appendix D: Einstein crystal 


For thermodynamic integration we are interested in the 
limit A: —>■ 0. In this limit there is no penalty for moving 
the system as whole, hence the mean squared displace¬ 
ment becomes of the order of , where L is the box side 
length. This result means that the function (jr — rgj^)fc 
will be strongly peaked at fc = 0, thus making its integra¬ 
tion difficult. For this reason, we would like this function 
to vary slowly with k. This behaviour can be achieved by 
fixing the centre of mass of the system, so that drifting 
as a whole is forbidden [107]. 

The centre of mass is defined as: 


rcM = 


where fXi 


rui 


(D5) 


When computing the potential energy for the harmonic 
spring, we must apply the following correction: 

N 


| r ( C )_ r (^)|2 




fCM)|2 


(D6) 


I 


The harmonic potential is defined as follows: 

y(r|rg,A:) = ^[r-rgj^ = (Dl) 

i 


where i is the index for the i-th. particle and C and 
U denote the corrected and the uncorrected coordinates 
respectively. The configurational partition function re¬ 
quires a correction, hence we define the corrected parti¬ 
tion function Z^yi with centre of mass fixed at rcM = 0 
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and note that: 



(D7) 


where solution of the integral was obtained after a fair 
amount of algebra by rewriting the Dirac delta as the 
Fourier sum d(x) = l/(27r^) f dkexp(ikx) [103, 118]. 

Using Eq. (D4) we find the mean squared displacement 
for the constrained Harmonic oscillator: 


„ „ ,2^ {N-l)d 

r-ro|)cM = 2 -—- =-—-. 


\ dk J 


NVT 


pk 


(D8) 


This result can be interpreted as the mean squared dis¬ 
placement of the (TV — l)d harmonic oscillator: fixing 
the centre of mass is equivalent to fixing one particle 
and integrating Eq. (D7) over the remaining degrees of 
freedom by doing the change of variables r' = — rjv 

(conveniently with unit Jacobian) if the iV-th particle is 
fixed. 

To conclude, let us relabel the potential as 
U(r|ro,A:, A) = (1 - A)$(r) -h ^Xk\r-rof, (D9) 


where 4>(r) is an arbitrary field, it could be, for instance, 
an additional inter-atomic interaction independent of k 
or, even the zero field. Let us consider the limit A —>■ 0: 
from the ratio of the partition functions for the con¬ 
strained and unconstrained centre of mass, we find: 


Appendix E: Polydisperse hard-sphere fluid and 
total accessible volume 

We can write the total accessible volume as 


- \nV{N,p) = -iVlnUbox + A'/ex(</'), (El) 


where (j) is the volume fraction and fex( 4 >) is the excess 
free energy, which is the difference in free energy between 
the hard sphere fluid and the ideal gas. We can compute 
the excess free energy by thermodynamic integration 
[107]. We start by noting that 9F/cl(l/Ubox) = 
and define the number density p = iV/Vbox, hence we 
write 


/ex(p) 


F{p) 

N 


E(‘d)(p) 

N 





(E2) 

By noting that the volume fraction of a polydisperse sys¬ 
tem \s (p = Vdp{<T'^) [119], where Vd is the volume of the 
d-dimensional unit sphere and {a‘^) is the d-th moment 
of the distribution of diameters, we can change variable 
and write 


/ex(<(>) = 


F{P) 

N 


N 





(E3) 

where = P/p is the compressibility factor (we set 
/3 = 1 everywhere). 

Analytical approximations for the compressibility fac¬ 
tors for the two and three-dimensional polydisperse hard 
sphere fluid have been proposed. For the hard disk fluid 
we use the Santos-Yuste-Haro (eSYH) equation of state 
[119] 


Zcm{X = 0 ) 

Z{\ = 0) 



(DIO) 



'U(rcM = 0), 


where <5 is the Dirac delta function and V{ycm — 0) is 
the probability density of the centre of mass being at 0 
when A = 0. Hence we write: 

^CM(A = 0) = Z(A = 0)iP(rcM = 0) (Dll) 
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(E4) 


where c/q = -k/ y/V2 is the crystalline close packing frac¬ 
tion. 

For three-dimensional fluids, depending on the volume 
fraction, we choose two different equations of state. For 
volume fraction 0 > 0.5, Santos et al. [119] suggest 
the following equation of state based on the Carnahan- 
Startling (CS) equation of state for the monodisperse 
fluid: 


where V depends on the details of the system. If the 
equilibrium structure is invariant to translations, a con¬ 
dition that holds true in a system with periodic boundary 
conditions, then we can take P = l/14eii, where Well is 
the smallest repeating unit in the periodic system (unit 
cell). This is at worst Well = Wox, while for a fee Ein¬ 
stein crystal it would correspond to the Wigner-Seitz cell 
Well = Wox/TV [102]. 


= 1+ 
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For volume fractions 4> < 0.5 the eCSK equation 
of state should be preferred (based on the Carnahan- 
Starling-Kolafa equation of state for the monodisperse 
fluid) 


^ 6(a3)2 ))■ 


The excess free energy can thus be obtained by substi¬ 
tuting one of Eq. (E4) to (E6) in the integral of Eq. (E3) , 
which can then be evaluated numerically for the desired 
volume fraction. 


which is a unimodal distribution with mode Xm = 
exp(A^cr^/K -I- A* — cr^). The distribution is such that 


POO 

/ I{x] fi,a, N) dx = exp 
Jo 


'a^N^ 

2k2 


Nn 

At 


(F2) 


Since a oc l/N, for large N we have I{x\ p,, a, N ^ 
At) dx = £exp(7V/t/At), with s some constant. Thus in 
the thermodynamic limit {N, V,lja ^ oo) we obtain the 
expression for the Gibbs configurational entropy per par¬ 
ticle, see Eq. (12). 


Consider Eq. (6) which we rewrite in terms of the di- 
Appendix F: Thermodynamic limit of the mensionless free energy per particle 

generalised configurational entropy 


Erom the fit to the empirical c.d.f. of B{'P) with the 
generalised log-normal cumulative distribution function, 
corresponding to Eq. (28), we obtain the set of param¬ 
eters /r, cr, and C- From the inset in Fig. 5, we observe 
that the mean p, and scale parameter cr scale linearly with 
\/N. In particular we note that a seems to approach 
zero in the thermodynamic limit, as expected. Further¬ 
more we note that the shape parameter C, seems to be 
approximately independent of l/N and to have a value 
of approximately 2 for all system sizes, thus suggesting 
that the distributions of pressures are consistent with a 
log-normal distribution. 

Therefore, under the reasonable assumption that the 
biased distribution of pressures B{x\V) is log-normal, we 
write the integrand in Eq. (11) as 


I{x;fi,a,N)= B{x\V)x^/'^ = 

1 


!;V2 


: exp 


TTa^ 




(Fl) 


/(iP|</'Ss) =-^ln('y) = c-H-ln('P), (F3) 

iV AC 


where c = C{N)/N; to test the consistency of the results 
thus obtained we compare the expression for Edwards 
configurational entropy, see Eq. (12), to the Gibbs con¬ 
figurational entropy per particle, see Eq. (23). We thus 
find that 


(/)b =c + -(ln('P))B, 

K, 


(F4) 


which is consistent with Eq. (E3) in the thermodynamic 
limit. We have thus correctly recovered the power law 
relation between pressure and basin volume, Eq. (6). 
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